
% 模拟退火算法求解非线性规划问题
% 目标函数: f(x,y) = sin(2x) + cos(2y) + ((x-1)^2 + (y-1)^2)/4
% 约束条件: 
%   -2 ≤ x ≤ 2
%   -2 ≤ y ≤ 2
%   x^2 + y^2 ≤ 3

function simulated_annealing_optimization()
    % 初始化参数
    rng('shuffle'); % 随机数种子，确保每次运行结果不同
    T_init = 100;   % 初始温度
    T_min = 1e-6;   % 终止温度
    alpha = 0.95;   % 温度衰减系数
    max_iter = 100; % 每个温度下的迭代次数
    
    % 定义约束边界
    x_min = -2; x_max = 2;
    y_min = -2; y_max = 2;
    
    % 生成初始解 (在可行域内随机生成)
    current_sol = generate_initial_solution(x_min, x_max, y_min, y_max);
    current_val = objective_func(current_sol(1), current_sol(2));
    
    % 记录最优解
    best_sol = current_sol;
    best_val = current_val;
    
    % 模拟退火主循环
    T = T_init;
    while T > T_min
        for i = 1:max_iter
            % 生成新解 (在邻域内随机扰动)
            new_sol = generate_neighbor(current_sol, T, x_min, x_max, y_min, y_max);
            new_val = objective_func(new_sol(1), new_sol(2));
            
            % 计算能量差
            delta_E = new_val - current_val;
            
            % 决定是否接受新解
            if delta_E < 0
                % 新解更好，直接接受
                current_sol = new_sol;
                current_val = new_val;
                
                % 更新全局最优解
                if new_val < best_val
                    best_sol = new_sol;
                    best_val = new_val;
                end
            else
                % 新解更差，以一定概率接受(避免陷入局部最优)
                p = exp(-delta_E / T);
                if rand() < p
                    current_sol = new_sol;
                    current_val = new_val;
                end
            end
        end
        
        % 降低温度
        T = T * alpha;
        
        % 显示当前信息
        fprintf('温度: %.4f, 当前解: (%.4f, %.4f), 当前值: %.4f, 最优值: %.4f\n', ...
                T, current_sol(1), current_sol(2), current_val, best_val);
    end
    
    % 输出最终结果
    fprintf('\n优化结果:\n');
    fprintf('全局最优解: x = %.6f, y = %.6f\n', best_sol(1), best_sol(2));
    fprintf('最优目标函数值: %.6f\n', best_val);
    
    % 绘制目标函数曲面和最优解
    plot_objective_function(best_sol);
end

% 目标函数定义
function f = objective_func(x, y)
    f = sin(2*x) + cos(2*y) + ((x-1)^2 + (y-1)^2)/4;
end

% 生成初始解 (确保在可行域内)
function sol = generate_initial_solution(x_min, x_max, y_min, y_max)
    while true
        x = x_min + (x_max - x_min) * rand();
        y = y_min + (y_max - y_min) * rand();
        if x^2 + y^2 <= 3  % 检查是否满足约束条件
            sol = [x, y];
            return;
        end
    end
end

% 生成邻域解
function new_sol = generate_neighbor(current_sol, T, x_min, x_max, y_min, y_max)
    % 温度越高，扰动范围越大
    scale = T / 100;
    
    while true
        % 在当前解附近随机扰动
        x = current_sol(1) + scale * (2*rand()-1);
        y = current_sol(2) + scale * (2*rand()-1);
        
        % 确保新解在可行域内
        x = max(x_min, min(x_max, x));
        y = max(y_min, min(y_max, y));
        
        if x^2 + y^2 <= 3  % 检查是否满足约束条件
            new_sol = [x, y];
            return;
        end
    end
end

% 绘制目标函数曲面和最优解
function plot_objective_function(best_sol)
    % 创建网格数据
    [X, Y] = meshgrid(linspace(-2, 2, 50), linspace(-2, 2, 50));
    
    % 计算目标函数值
    Z = sin(2*X) + cos(2*Y) + ((X-1).^2 + (Y-1).^2)/4;
    
    % 绘制曲面
    figure;
    surf(X, Y, Z);
    hold on;
    
    % 标记最优解
    best_z = objective_func(best_sol(1), best_sol(2));
    plot3(best_sol(1), best_sol(2), best_z, 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r');
    
    % 图形设置
    xlabel('x');
    ylabel('y');
    zlabel('f(x,y)');
    title('目标函数曲面及最优解');
    colorbar;
    shading interp;
    view(-45, 30);
    
    % 绘制约束边界
    theta = linspace(0, 2*pi, 100);
    r = sqrt(3);
    x_circle = r * cos(theta);
    y_circle = r * sin(theta);
    z_circle = sin(2*x_circle) + cos(2*y_circle) + ((x_circle-1).^2 + (y_circle-1).^2)/4;
    plot3(x_circle, y_circle, z_circle, 'k-', 'LineWidth', 2);
    
    hold off;
end